%%% Written by  Daniel Lewis, 2020, for "Robust Inference in Models 
%%% Identified via Heteroskedasticity", Review of Economics and Statistics.
%
% For a 2-variable model, this function takes as inputs the Tr by 3
% matrices of the unique vectorization of outer products of residuals for 
% each regime, the total sample length, the nuisance parameters, the null 
% for the parameter of interest, and the index of the parameter of 
% interest. It returns the GMM objective function at the supplied 
% parameters using the efficient weighting matrix, the sum of moments, the 
% covariance of moments, and the efficient weighting matrix. This is for 
% use in computing projection or subset tests following Theorem 3, where 
% Hint is the parameter of interest.
%
% Inputs:
% etasq1:   outer product of regime 1 residuals (T1x3 unique vectorization)
% etasq2:   outer product of regime 2 residuals (T2x3 unique vectorization)
% T:        total sample length
% nuis:        vector of nuisance parameters
% Hint:     null value of parameter of interest, H_21 or H_12
% pind:     index of parameter of interest, 1 for H_21, 2 for H_12
% Outputs:
% obj:      GMM objective function
% f_T:      sum of moments
% matinv:   efficient weighting matrix
% omega:    moment covariance matrix
function [obj,f_T,matinv,omega]=momfastsub(etasq1,etasq2,T,nuis,Hint,pind)
p=NaN(6,1); el=1:6; el(pind)=[]; p(pind)=Hint; p(el)=nuis; 
V=[1,p(2);p(1),1]*diag([p(3),p(4)])*[1,p(2);p(1),1]'; % construct fitted reduced form covariance matrix for regime 1
V1=[V(1:2),V(4)]; % vectorize unique elements
V=[1,p(2);p(1),1]*diag([p(5),p(6)])*[1,p(2);p(1),1]'; % construct fitted reduced form covariance matrix for regime 2
V2=[V(1:2),V(4)]; % vectorize unique elements
f_T=[[etasq1-V1;zeros(length(etasq2),3)] [zeros(length(etasq1),3); etasq2-V2]]; % compute moments and stack
omega=f_T'*f_T/T; % compute covariance matrix of moment equations
f_T=sum(f_T)'; % compute sum of moments
matinv=nearestSPD(omega^-1); % compute efficient weighting matrix (nearestSPD corrects in case of numerical issues with symmetry or positive definiteness)
obj=f_T'*matinv*f_T/T;
% If some numerical error results in a negative objective function in spite
% of nearestSPD, set objective function to a large value to continue search
% in another direction.
if obj<0
    obj=1e8;   
end
end
